#include <stdio.h>
#include <cuda.h>
#include <sys/time.h>
#include <cuda_runtime.h>
#include "cublas_v2.h"

#define M 2048
#define N 2048
#define K 2048

# define timersub(a, b, result)                                         \
	do {                                                            \
		(result)->tv_sec = (a)->tv_sec - (b)->tv_sec;           \
		(result)->tv_usec = (a)->tv_usec - (b)->tv_usec;        \
		if ((result)->tv_usec < 0) {                            \
			--(result)->tv_sec;                             \
			(result)->tv_usec += 1000000;                   \
		}                                                       \
	} while (0)

__global__ void kernel2(float *C, const float *A, const float *B, int m, int n,  int k)
{
	const int row = blockIdx.y * blockDim.y + threadIdx.y;
	const int col = blockIdx.x * blockDim.x + threadIdx.x;

	if (row >= m) return;
	if (col >= n) return;

	float sum = 0;
	for (int i = 0; i < k; ++i) {
		sum += A[row * k + i] * B[i * n + col];
	}

	C[row * n + col] = sum;
}

int main(int argc, const char *argv[])
{
	struct timeval time_start, time_end, time_res;

	int count = 0;
	float alpha = 1.0;
	float beta = 0;

	float *m1_h, *m2_h, *o_h;
	float *m1_d, *m2_d, *o_d;
	float *o_h2;

	m1_h = (float *) malloc(M * K * sizeof(float));
	m2_h = (float *) malloc(K * N * sizeof(float));
	o_h = (float *) malloc(M * N * sizeof(float));
	o_h2 = (float *) malloc(M * N * sizeof(float));

	cudaMalloc( (void **)&m1_d, M * K * sizeof(float));
	cudaMalloc( (void **)&m2_d, K * N * sizeof(float));
	cudaMalloc( (void **)&o_d, M * N * sizeof(float));

	// Initialize host array and copy it to CUDA device
	for (int r = 0 ; r < M; ++r) {
		for (int c = 0 ; c < K; ++c) {
			m1_h[r*K + c] = rand()/(float)RAND_MAX;
		}
	}

	for (int r = 0 ; r < K; ++r) {
		for (int c = 0 ; c < N; ++c) {
			m2_h[r*N + c] = rand()/(float)RAND_MAX;
		}
	}

	cudaMemcpy(m1_d, m1_h, M * K * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(m2_d, m2_h, K * N * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(o_d, o_h, M * N * sizeof(float), cudaMemcpyHostToDevice);

	cublasHandle_t handle;

	if(cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS){
		printf("cublasCteate fail\n");
		exit(-1);
	}
	for(int i = 1 ; i <= 16 ; i++){
	//set grid size and block size
//	int i = 16;
	dim3 grid((M + i - 1)/i, (N + i - 1)/i, 1);
	dim3 threads(i, i, 1);

	//timer start
	gettimeofday(&time_start, NULL);

	kernel2 <<< grid, threads >>> ( o_d, m1_d, m2_d, M, N, K );

	cudaDeviceSynchronize();

	//timer end
	gettimeofday(&time_end, NULL);

	timersub(&time_end, &time_start, &time_res);

	printf("%ld,%ld ms\n", time_res.tv_sec, time_res.tv_usec/1000);

	}
	cudaMemcpy(o_h,o_d,M*N*sizeof(float), cudaMemcpyDeviceToHost);
	cudaDeviceSynchronize();

	//get cublas result for comparison
	gettimeofday(&time_start, NULL);

	cublasSgemm(handle, CUBLAS_OP_T,CUBLAS_OP_T, M, K, N, &alpha, m1_d, K, m2_d, N, &beta, o_d, N);
	cudaDeviceSynchronize();

	cudaMemcpy(o_h2, o_d, M*N*sizeof(float), cudaMemcpyDeviceToHost);

	//compare differences (transpose for colomn major)
	for(int i = 0 ; i < M ; i ++){
		for(int j = 0 ; j < N ; j++){
			if( (o_h[i*N+j] - o_h2[j*N+i]) > 0.001 || (o_h[i*N+j] - o_h2[j*N+i]) < -0.001){
				count++;
			}
		}
	}
	printf("number of differences : %d\n",count);

	free(m1_h);
	free(m2_h);
	free(o_h);

	cudaFree(m1_d);
	cudaFree(m2_d);
	cudaFree(o_d);
}

